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Abstract: In this article we describe simulations of Higgs boson production via the gluon 
fusion and Higgs-strahlung processes, using the positive weight next-to-leading-order (NLO) 
matching scheme, POWHEG, in the Herwig++ 2.3 event generator. This formalism consis- 
tently incorporates the full NLO corrections to these processes within the parton shower 
simulation, without the production of negative weight events. These simulations include a 
full implementation of the truncated shower required to correctly model soft emissions in 
an angular- ordered parton shower. We present a thorough validation of these simulations, 
comparing them to other methods and calculations. The results obtained for the gluon 
fusion process corroborate, and provide detailed explanations for, findings reported by Ali- 
oli et al using an independent POWHEG simulation, neglecting truncated shower effects, 
released at the same time as our code. 
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1. Introduction 

The primary objective of the ongoing and imminent physics programmes at the Tevatron 
and LHC, is to elucidate the nature of electroweak symmetry breaking. The great majority 
of the effort in this direction is devoted to the hunt for the Higgs boson, the origin of this 
symmetry breaking in the Standard Model (SM) [1-4]. 

Of all the ways in which the SM Higgs boson can be produced, the gluon fusion process 
[5], in which it couples to colliding gluons via a top quark loop, has the largest cross section 
for Higgs boson masses less than ~1 TeV. This process is of great importance for the 
detection of the Higgs boson at the Tevatron and, more so, at the LHC, particularly in the 
low mass region, favoured by the latest fits of the Standard Model to electroweak precision 
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data [6], where the decay of the Higgs boson into two photons is known to give a clean 
experimental signal. Although observing a narrow resonance in the diphoton invariant mass 
spectrum should be possible using only the experimental data [7], determining the quantum 
numbers and couplings of the resonance i.e. determining that it really is a fundamental 
scalar and, moreover, the SM Higgs boson, will involve a comprehensive analysis of a number 
of channels, using accurate, flexible, Monte Carlo simulations to predict distributions for 
signals and backgrounds. 

A key part of that identification procedure will be the measurement the couplings of the 
Higgs boson to the weak gauge bosons and the top quark. Although the gluon fusion process 
directly probes the latter, associated QCD radiation renders it a significant background to 
the vector boson fusion process (VBF), in which the Higgs boson originates from HWW 
and HZZ vertices. Precise simulations of the gluon fusion process are then also required 
to model the extent to which these events contaminate the VBF signal [8-10]. 

Another direct probe of Higgs boson interactions with electroweak gauge bosons is the 
Higgs-strahlung process [11,12]. At leading order this consists of a quark anti-quark annihi- 
lation producing a virtual vector boson (V = W/Z), which becomes on-shell by radiating a 
Higgs boson. This channel is particularly important for Higgs boson searches at the Teva- 
tron [13]. The cross section for Higgs-strahlung processes at the LHC is approximately one 
order of magnitude larger than at the Tevatron, however, the backgrounds scale up by a 
much greater factor. In spite of these difficulties the Higgs-strahlung processes can still be 
observed, notably in the 77 and ly+iy - decay channels, where significant background re- 
jection can be achieved (see for example Ref. [14] for recent experimental studies). Interest 
in the bb plus leptons decay channel has also recently been revived through the observa- 
tion that a highly collimated bb jet, from a boosted Higgs boson, could provide a clean 
signature [15,16]. These processes were recommended for inclusion in studies to determine 
Higgs boson couplings in Ref. [17]. As with the gluon fusion process, accurate Monte Carlo 
simulations will be central to such studies. 

Lately combined Tevatron analyses of the gluon fusion and Higgs-strahlung channels 
have begun to show sensitivity to Higgs boson production [18]. Presently these studies 
exclude, at the 95% confidence level, the existence of a SM Higgs boson with a mass in 
the range 160 - 170 GeV. This exclusion limit largely follows from analysis of the gluon 
fusion channel in which the Higgs decays into a W + W~ pair. Monte Carlo simulations of 
Higgs-strahlung and gluon fusion processes were essential in obtaining these measurements. 

In recent years research in Monte Carlo simulations has seen major progress, most 
significantly in the extension of existing parton shower simulations to consistently include 
exact next-to-leading order corrections [19-34] and, separately, in the consistent combi- 
nation of parton shower simulations and high multiplicity tree-level matrix element gen- 
erators [35-40]. The first successful NLO matching scheme was MC@NLO [19-24] which 
has been realised using the HERWIG event generator for many processes. This method 
has two draw backs; first, it involves subtracting the parton shower approximation from 
the NLO calculation, which can lead to unphysical negative weight events, second, the 
implementation of the method is fundamentally dependent on the details of the parton 
shower algorithm. In 2004 a new formalism known as POWHEG (POsitive Weight Hardest 
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Emission Generator [26]) was derived, achieving the same aims as MCONLO but with the 
added benefits of only generating physical, positive weight, events as they occur in na- 
ture, and of being independent of the details of the parton shower algorithm. This method 
has been already successfully applied to a number of phenomenologically important pro- 
cesses [25,27,30-34,41,42]. 

In this paper we describe the application of the POWHEG method to Monte Carlo sim- 
ulations for Higgs boson production via gluon fusion and Higgs-strahlung processes, within 
the Herwig++ [43, 44] event generator. We include a complete description of truncated 
shower effects. 1 Our primary aim is to present the ingredients used in the simulations and 
to validate them, where appropriate, against existing calculations. 

The structure of the paper is as follows. In Sect. [2] we briefly review the main fea- 
tures of the POWHEG method. In Sect. || we collect the essential formulae relating to the 
NLO cross sections, for implementation in the program. In Section |] we give details of the 
event generation process for the hard configurations. This is accompanied by a descrip- 
tion of how the hard configurations are subsequently reproduced by the angular- ordered 
parton shower, including truncated showering effects, in respect of the colour coherence 
phenomenon manifest in soft wide angle gluon emissions. In Sect. |5| we present the re- 
sults of our implementation, comparing it to the MCFM and MCONLO generators, before 
summarizing our findings in Sect. [f| 



2. The POWHEG formalism 



The central formula in the POWHEG approach is derived by manipulating and modifying, 
through the inclusion of formally higher order terms, the NLO differential cross section, 
such that it has the same form as that given by the parton shower [26]. For a given iV-body 
process this is 



dcr = B d$_e 
where B (<3?_b) is defined as 

B($ B ) = B{<& b ) + V{<S>b)+ I d$ 



R(®R 

A « (0) + m r A * {kT (<&B ' d * R 



(2.1) 



R($B,®R)-J2 C i(® B >®^ 



(2.2) 



with B(<&b) the leading-order contribution, dependent on the Born variables, <3?#, which 
specify a point in the iV-body phase space. V(&b) is a finite contribution arising from 
the combination of unresolvable, real emission and virtual loop corrections. The remaining 
terms are due to the iV+l-body real emission processes, hence they have an additional 
dependence on the radiative variables, $r, which parametrize the phase space associated 
with the extra parton. The real emission term, R(&b, &r), is given by the product of the 
parton flux factors with the relevant squared real emission matrix element, summed over 



1 Currently, only the Herwig+- \- POWHEG simulation of the Drell-Yan process includes a complete treat- 
ment of truncated shower effects [32]. 
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each channel contributing to the NLO cross section. C% (&b,^r) denotes a combination of 
real counterterms / counter-event weights, regulating the singularities in R (<&b, &r)- Finally, 
the POWHEG Sudakov form factor, A^, is defined as' 



A i? (Pt) = exp 



(2.3) 



where kx (&b, $r) tends to the transverse momentum of the emitted parton in the soft and 
collinear limits. 



To O (as) Eq.^jJ is just the usual NLO differential cross section. The analogous 
parton shower expression is given by replacing B(&b) - > B($b) and the real emission 
corrections R(&Bi&r) by their collinear approximation, i.e. replacing the argument of 
the Sudakov form factor by a sum of Altarelli-Parisi kernels. Typically parton shower 
simulations generate an iV-body configuration according to B ($_b) and then shower it using 
such a Sudakov form factor. In the POWHEG formalism the initial iV-body configuration 
is instead generated according to B(&b), and retained with probability (0) as a non- 
radiative event, or, showered to give the hardest emission with px = kx ($b, $r), with 
probability A r (pt) (Eq. |2.ip . Further, lower pr, emissions represent terms of next-to- 
next-to-leading-order (NNLO) and beyond, hence we can return to the usual parton shower 
formulation to simulate these. Provided the perturbative approach is valid, i.e. provided 
the next-to-leading order terms are smaller than the leading order ones, it is clear from 



Eq. 2^2 that B ($_b) is positive, therefore no negative weights arise. 

Furthermore, it is well known that when a bunch of collimated QCD charges emit 
a gluon at wide angle, the intensity of the radiation is proportional to the coherent sum 
of emissions from the constituents, viz. the jet parent. This effect is manifest in the 
perturbative series as large soft logarithms. The other major triumph of the POWHEG [26] 
approach, besides avoiding negative weights, is in rigorously decomposing the angular- 
ordered parton shower into a truncated shower, describing soft wide angle radiation, the 
hardest emission, as described above, and a further vetoed shower comprised of lower pt, 
increasingly collinear emissions. 

We implement the POWHEG formalism in fullness according to the following procedure: 



generate N- and iV+l-body configurations according to Eq. 2.1 
directly hadronize any iV-body, non-radiative, events; 

map the radiative variables parametrizing the emission into the evolution scale, mo- 
mentum fraction and azimuthal angle (qh, Zh, fa), from which the parton shower 
would reconstruct identical momenta; 

take the initial iV-body configuration, generated from B($>b), and evolve the 
emitting leg from the default initial scale down to using the truncated shower; 

insert a branching with parameters (qh, z^, fa) into the shower when the evolution 
scale reaches q^; 

generate pr vetoed showers from all external legs. 
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3. Next-to-Leading Order Cross Sections 

The NLO ingredients needed to implement the POWHEG method for Higgs-strahlung pro- 
cesses can be performed in the same way as Ref. [32] due to the careful organization of the 
Drell-Yan differential cross section in our earlier work. There we explicitly factorized the 
real emission corrections to the leading-order q + q — ► I + 1 process at the level of the phase 
space and the real emission matrix elements. This meant that the radiative variables could 
be generated completely independently from the details of the decay of the off-shell vector 
boson. We refer the reader to Ref. [32] for details of those matrix elements. 

Since the diagrams involved in the NLO corrections to Drell-Yan and Higgs-strahlung 
processes are identical up to replacing the final-state lepton pair with a vector boson and 
a Higgs boson, the factorized NLO differential cross section is exactly the same as that in 
Ref. [32] for the Drell-Yan process; one simply replaces the q + q — > I + I leading-order 
matrix element with that for q + q — > V + H throughout. No information regarding spin 
correlations is lost in this procedure and the full NLO distribution is generated without 
approximation. This method is originally due to Kleiss [45,46] and was extended for use in 
our POWHEG simulation, with a view to making other processes easier to implement. Due 
to the complexity of the 4-body final state (we include the decays of V and H) and the 
subtleties of the Monte Carlo algorithm, we have carried out detailed comparisons against 
the parton level NLO simulation MCFM [47]. Results of these comparisons are given later 
in Sect.g. 

The NLO cross section for Higgs boson production via gluon fusion was first calculated 
in Refs. [48,49], in the infinite top-quark mass limit. Later this calculation was improved 
to exactly include the effects of the finite top mass in [50], where it was found that the 
infinite top-quark mass limit provided an excellent approximation to the full result; below 
the ti threshold the exact and approximate calculations agree at the level of 10% [7]. In the 
last decade the NNLO corrections were also fully computed by four independent groups, 
in the infinite top-quark mass limit [51-55]. More recently a parton level Monte Carlo 
program accurate to NNLO has become available [56,57]. A key point arising from this 
theoretical activity is that the dominant perturbative corrections to the total cross section 
originate from virtual and soft-gluon corrections, which explains the surprising accuracy of 
the infinite top-quark mass approximation for ma < 2m t [58]. 

In our simulation we use the infinite top-quark mass limit. Although the NLO formulae 
in [48,49] are helpful, they were only used to compute the most inclusive of measurements, 
namely the total cross section. Further work is needed to obtain a form suitable for a fully 
differential Monte Carlo simulation. In this section we collect the ingredients that arise in 
the NLO calculation for g + g — > H, necessary for the implementation of the POWHEG 
method. 

3.1 Kinematics and phase space 

Here we restrict ourselves to considering leading-order processes of the type, 
Pe> + Po —* Pi + ••• +PN, hi which all the particles in the iV-body final state are either 
massive or colourless. For such processes the NLO corrections may contain soft singulari- 
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ties and initial-state collinear singularities only. The incoming hadron momenta are labeled 
P@, for hadrons incident from the ±z directions, respectively. It therefore follows that the 
momenta of the colliding, massless partons, with momentum fractions and Xq, are given 
by p @ = x @ -P @ . The momentum of the ith final-state parton produced in the leading-order 
process is denoted pi. 

Using p to represent the sum of all pi, and y to signify the rapidity of p, the phase 
space for the leading-order process can be written simply as 

d&B = dx©dx e dl>B = -dp 2 dydl>B, (3.1) 

where d<3>_B is the Lorentz invariant phase space for the partonic 2 — > TV process and s the 
hadronic centre-of-mass energy. We refer to the variables = |p 2 5 y, as the Born 

variables, where <&b parametrizes the iV-body phase space in the partonic centre-of-mass 
frame; for the simple case of a decaying scalar, like the Higgs boson, the matrix element 
does not depend on <3?#, which could therefore be trivially integrated out at this point. 

The NLO real emission corrections to the leading-order process consist of 2 — > N + 1 
processes, p®+Pe — > pi + ...+pN + k, where we denote the momenta of the same N particles 
produced in the leading-order process pi and that of the extra colour charged parton by k. 
For these processes we introduce the Mandelstam variables s, t, u and the related radiative 
variables &r = {x, y}, which parametrize the extra emission: 

s = (Pe +Pef = — , (3-2a) 
x 

i=( P(B -k) 2 =- lP l(l- X )(l-y), (3.2b) 

u=(p e -k) 2 =-^(1-i)(1 + j,), (3.2c) 



where p = Yl Pi- We do not explicitly include an azimuthal angle for the gluon about the 
z axis, instead it is used to define the +y axis relative to which the azimuthal angle of the 
other final-state particles is measured; ultimately all generated events are randomly rotated 
about the z axis in the hadronic centre-of-mass frame. 

To perform a simultaneous Monte Carlo sampling of the N- and iV+l-body phase 
spaces one has to specify the integration variables. We choose two of these to be the mass 
and rapidity of the system of colourless particles, therefore p 2 = p 2 and y = y, where y 
is the rapidity of p. 2 An immediate consequence of this partial mapping of the N and 
iV+1 body phase spaces is that the momentum fractions, rc@, of the incident partons in the 
2 — > N + 1 processes are related to those of the 2 — > N process by 



*e / 2-(l-x)(l-y) _xq \ 2 - (1 - x) (1 + y) 

e " v^p-(i-x)(i + y)' Xe ~ ^p-(i-x)(i- y y 

where, by definition, p@ = x^P^. 



2 Henceforth we will always refer to these variables as p 2 and y. 
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Since we restrict ourselves to processes for which the NLO corrections contain at most 
soft and initial-state collinear singularities, the product of tu with the squared real emission 
matrix elements will be finite throughout the radiative phase space. Working in conventional 
dimensional regularization, in n = 4 — 2e space-time dimensions, we are then able to 
perform an expansion in e of the iV+l-body phase-space measure, d$Ar+i, using similar 
manipulations to those in Refs. [59-62], giving 

dSjy+i = d$ B d$ fl P (\) crJ(x,y), (3.4) 
(47r) x z \P J 

where here the Born variables <&b specify a configuration in the rest frame of p rather than 
p. The function J (x, y) is given by 

iu 

J(x,y) = [S5(l-x)+C(x) (26(l+y) + 26(l-y))+H(x,y)] ^ , 

where 

S = \ --ln7] + 8ln 2 r], (3.5a) 



Cx = — - --- -lnx + 2 — '- , 3.5b 

e(l-aO p 0--x) \ 1-x J 



H{x,y) = 



(1-x), 



1 \ ( 1 
+ 



1 + yJ+ \ 1 -y/ + 



(3.5c) 



with p = (^2i yP^^j /$ an( l V = v 7 ! — P- The constant cp, which appears due to the use 
of dimensional regularization, is given by 

cr = U.Y r(1 " e)2r(1 + £) 
r 1 > r(l-2e) 

and the ^-distributions are defined by the relation 

nl / ln n (l-x)\ /- 1 ln n (l-x) 



dx/t(x) ( ln "/ 1 _ a; X) ) =/ da 



p \ - * / p j p I x 

for any sufficiently regular test function h (x) . The radiative phase space can be parametrized 
in terms of the radiative variables as 

d$ K = idydx. (3.6) 

The precise definition of $s in the context of the radiative event is, in general, given by 
a series of boosts and rotations, denoted by B, which embed the iV-particle, leading-order 
final state in the iV+l-particle radiative events. To assemble a radiative event configuration 
we first construct the N final-state momenta of the leading-order configuration according to 
the definition of &b and, separately, the momenta of the incident partons and the radiated 
parton k in the hadronic centre-of-mass frame: 

P® = ^Vs{xe,0,0,+x & ) , p= (£ T coshy, 0, -p T , ^sinhy) , (3.7) 

Pe = 7^(^9,0,0, -x e ) , k = (p T coshy fc , 0, p T , Pt sinhy fe ) , (3.8) 
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where Et = yPx +P , and 

2 P 2 l h 2Wl ^2 , 1, f 1 + v\^ 1 , ( 2-{l-x)(l-y) 

(3.9) 

To embed the iV final-state particles of the leading-order process in the radiative event we 
first transform them with a longitudinal boost, B y , taking us to their rest frame. Should 
we wish to adhere to the conventions in Refs. [59-62], we then apply a rotation to them, R, 
defined such that, ultimately, B will preserve the direction of p®, i.e. with p e defining the 
+z axis in the p rest frame, and also with the transverse momentum of k defining the y axis. 
However, for the simple case at hand, no such rotation is necessary since the Higgs boson 
decays isotropically in its rest frame. 3 A further transverse boost, By, parallel with the y 
axis, is then carried out, where B^ 1 is a transformation from the lab to the frame in which 
p has no transverse component. Finally, the inverse boost B~ 1 is applied to the N particles, 
which clearly returns them to a frame where their total rapidity is y = y. Altogether the 
embedding boost B, to be applied to the leading-order final-state momenta, is given as 

B = By 1 B T MB y , (3.10) 



which, combined with p®, p Q and k in Eqs. |3/7|, |3,8|, completely specifies the radiative 
kinematics. 

3.2 Differential cross section for a + b — ► n 

In this section we restrict ourselves to discussing the general form of the NLO differential 
cross section for processes of the type a + b — ► i\ + ... + ijv, where n = Y^j ij are colourless 



particles, which we collectively refer to as neutrals. In Sect.^OI we give the squared matrix 
elements for g + g — > H which were ultimately inserted in these formulae. We reiterate 
that analogous Higgs-strahlung ingredients are identical to those in Ref. [32] modulo the 
substitution of the q + q — ► I + I leading-order matrix element with that of a given Higgs- 
strahlung process. 

We define the combined incident flux of partons of type a from hadron A, and partons 
of type b from hadron B, with respective momentum fractions x® and Xq, at a scale fx 2 , as 

£ab (x®,Xq) = fa {X(B,^ 2 ) fb { X Q, ^) , (3-11) 

where //(xj,^ 2 ) are the relevant parton distribution functions (PDFs). To obtain the 
differential cross section for a real emission process a + b — » n + c, we simply multiply the 



relevant squared matrix element by the phase-space measure in (Eq. 3.4), C a b and the flux 
factor l/2s : 

da - b+1 = 5§ M ab +1 (P®>Pe)£ab(x®,xe) d*W (3.12) 



3 As previously noted, the Higgs-strahlung process was simulated in a special way using the Kleiss trick 
as described in [32]. 
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In general each real emission process makes three contributions to the NLO differential 
cross section, one for each term in the phase-space measure Eq. 3J): a soft contribution 
der^ proportional to 5 (1 — x), collinear contributions dcr^ 0± proportional to S (1 ± y), and 
a finite, hard, contribution der^ proportional to Ti(x,y). The subscript in dcr^P and 
dcr^°® reflects the fact that they are bare, divergent, quantities. 

For leading-order processes of the type o + b — > n the two initial-state partons must be 
either a pair of gluons or a quark and an antiquark. The squared matrix elements for the 
real emission corrections to this process, in which a gluon is emitted from an initial-state 
parton, a + b — > n + g, factorize in the limit that the gluon is soft (x — ► 1), according to 

s2 



lim M^ b +1 (pe,p e ) = 8vrasM 



2c 



2C ab — M%(p e ,p e ) 



(3.13) 



The colour factor C a fc is equal to Ca if a and b are gluons, and Cjr if a is a quark and b is 
an antiquark or vice versa. The real emission processes a + 6 — > n + g are the only ones 
contributing to the cross section in the limit cc — ► 1, all other real emission matrix elements 
are finite in this limit, hence the product of tu with the squared matrix element vanishes 
there. Hence, the soft contribution to the NLO differential cross section is 



d*2 



a s c r ( H 



2ir 



p- 



7T 



^ °ab — --lnr/ + 161n 2 r, )B{$ B ) d<S> B 



(3.14) 



,e- 3 e 

where 5 ($b) is the differential cross section for the leading-order process 

d<Tab = B (* B ) , B (* B ) = ^2 (P©» Pe) ^afe (%, Xq) • 

For an arbitrary process in which an initial-state parton o, with momentum p, 
to produce a collinear time- like daughter parton with momentum k = (1 — £)p© and its 
space- like sister parton ac, with momentum xp@, the squared matrix element factorizes 
according to [63] 



(3.15) 
branches 



1 



s 2 



p 2 iu 



8iras/J> 



2c 



1 



P a& (x;e) M%(xp% 



Pe) 



(3.16) 



where we explicitly show the dependence of the iV and iV+1 parton squared matrix elements 
on the incident parton momenta. Replacing xp @ <-> p @ , a — > 6 in the splitting function we 
obtain the analogous formula for the case that c is emitted collinear to b (y — > — 1). Using 
this factorization of the matrix element, the collinear contributions to the real emission 
cross section, arising from initial-state radiation, are seen to have the general form, 



C% (x) 



+ ^ 



d^ @ , 



a s c r / n 



2tt 



P 



+ In 



P 



%c(P i i-c + 41n ^) 

£?6 x e) -d$eda:, 



B(* B ) d$ B , 



1 



(3.17) 



1 



■In 



P 

fJ> 2 X 



+ 2 



ln(l 



da 



CT® 
ab 



(1 - 

| g ^£ (*) s@ (*b) C.% (x e , are) ^ d<J> B dx, 



(1 - X) P.^ (x) - P £ r- (x), 
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where i = a in the case that parton a splits to produce parton c, and i = b for the case that 
parton b branches to produce c. P,r (x; e) and P.r- (x) are the bare and regularized Altarelli- 
Parisi kernels given in Appendix^; C-r- p.r- is equal to the coefficient of the 5 (1 — x) term 
in the latter, for the case p = 0. The B® and C® functions are related to the leading-order 
differential cross section and parton flux: 



B®(<5> B ) 
B e (<Z> B ) 



2p 
2p 



1 -M-ab ( X P®, Pe)£ab Offi, Xq) 



1 (p®, xp e ) C ab (x®, Xq) 



£-ab ( x ©' x e 



C ab (x 9 ,x e ) 



(3.18) 



IT 



In the MS scheme each collinear singular da ab term in the cross section is exactly 
compensated for by an additive collinear counter term identical to it (hence the CT label- 
ing). This amounts to absorbing the divergence in the PDFs, renormalizing them, hence 
we now omit them. The only remaining divergences are soft and collinear terms da ab , 
which we absorb in the soft contribution to the cross section d<7^. 

S C(±) 

Absorbing the soft and collinear terms da ab in the soft contribution to the cross 
section Eq. 3.14 , redefines it as 



d"2 



S B($ B ) d$ B 



(3.19) 



where 



So 



a s c r I ii 



2tt 



ab 



2 2 f P 2 \ , \ , 2 TT 2 ' 

^2 + -VaSg +ln I ~2 I (2p a Sg + 81nr?j + 16 In 77 - y 



(3.20) 



One must consider also the contribution of virtual corrections to the leading-order 
process as well as the real corrections above. Again, for a + b — ► n processes, these have a 
simple form (see e.g. Ref. [29]) 



x-*l 



Vo-M^Cp®, p, 



e) 



(3.21) 



with 



V 



2vr 



where .M^^ reg 



^2 ) ( '»'< 
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TT 

Y 



(p®, Pe) 



is the remainder of the virtual correction, regular as e 
leading-order squared matrix element, 



(3.22) 



0, divided by the 



M ah rcg (p e , Pe) - 



(p®, P 



'e 



Mf 6 (p®,Pe) 



(3.23) 



The poles in e in Vo will exactly cancel those in Sq making the differential cross section 
finite for e — ► 0. 

The cross section for the a + b — > n leading-order process, combined with the virtual 
corrections and a + b — ► n + g real emission corrections, may then be written as 



daab = B ($ B ) d$ B + V d$ B + -Rab (*B, d$B d$B, 



(3.24) 
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where V ($b) results from the combination of the soft and virtual corrections to the cross 
section (Eqs.|3,19|, |3720|) , 



V($ B ) = VB{<5> B ) 



V= 5o + Vo| 



e=0 ' 



(3.25) 



and R(&b,&r) results from the remaining collinear and hard real emission corrections, 



1 ^ 

Rab($B,$R) = —-n ab C ab { X(B ,XQ) B($ B ) , 
Z7T X 



T~tab 



2C®~ (x) 5(l-y) + 2C ft e - (x) 5 (1 + y) + H c 



(3.26) 



s 8vra 5 -M^(p©,p e ) 



H(x,y) . 



Note that in writing Eq. |3.24j we have tacitly equated M. (xp^,po) = M. (p®,xpQ 



M. (pq,pq) in the collinear term. This is possible due to our defining p = p and y = y 



in Section. |3.1| . 

The remaining contributions to the full NLO differential cross section are due to the 
production of n via new channels. For q + q — ► n processes we add contributions arising 
from the q + g — > n + q channel, 



flc 1 

R qg ($B,$R) = — -K qg £ qg (x®,X e ) , 



(3.27) 



n, 



H 



2C gq 6(l + y) + H qg , 
xiu 1 M^ g +1 (p®, p e ) 
s 8vtq 5 .M$,(p e , Pe) 



H(x,y) , 



and also a set of terms due to the q + g — > n + q channel, identical up to the replacement 
g -> g on all terms with the exception of C gg , For the 5 + g — > n process we add a 
contribution from g + q — ► n + q, for which has the same form as Eq. 3.27 , modulo 
interchanging the subscripts g<? <-> and a further contribution from the g + q ^ n + q 
process, derived from the former by replacing q — > q on all terms but C qg . 



3.3 g + g — » -ff matrix elements 

The squared, spin and colour averaged, leading-order matrix element for the g + g ^ H 
process is given by 

M» (P®, Pe) = M» = 576 ^_ e) , (3-28) 

a 2 n 2c 

where M = ^ a , with v the vacuum expectation value of the Higgs field and the scale 
// emerging from the use of conventional dimensional regularization. The real emission 
radiative corrections consist of three processes: g + g — ► H + g; q + g — ► H + q; and 
q + q ^ H + g. As noted above, the singular, soft and/or collinear, limits of the matrix 
elements, associated to the e poles in d$jv+ls are universal, therefore a full calculation of 
the matrix elements is only needed to multiply the TC (x, y) term in d$jv+i' Since the term 
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proportional to H(x,y) contains no poles in e, the full calculation of the squared matrix 
elements can be carried out in four dimensions, giving, 



M 



N+l 



3 1 



p 4 stu 



[ P s + s 4 + i A + u 4 



l <?9 



M N+l 
J l gq 



M N - +l 



(3.29) 



M N+l 
■ ' L gq 



The squared matrix element -M^ g +l can be obtained from -M^ q +1 by crossing symmetry 
simply by interchanging u t. 

The O (a$) virtual corrections to the g + g — > H process consist of purely gluonic 
swordfish and triangle vertex corrections, as well as a UV counter term. At NLO these 
corrections contribute to the cross section through their interference with the leading-order 
amplitude, their form is precisely that in Eq. 3.22 with, 



J l gg 



as 
2vr 



C A 
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4irb 
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99 



(3.30) 



where hr is the renormalization scale. When used with the general expressions in Sect.|3.2 



these are all the ingredients required to write down the B (<3?_b) function and the modified 
Sudakov form factors for this process (Eqs.^1, 2.3). 



4. Implementation 



4.1 Generation of the leading-order configuration 

The first phase of the simulation involves generating a leading-order configuration by sam- 
pling the B (&b) function (see Sect. ||, Eq. ^2), which is the NLO differential cross section 
integrated over the radiative variables, 



B(<S> B ) = B($ B ) 



ab 



Uj 1 ^ 

d^flTj K ab C ah (x®,xq) 

2tt x 



(4.1) 



Here the summation runs over all real emission processes a + b — > n + c contributing at 
NLO. Since the leading-order process is factorized inside the real emission terms R a b, the 
B (<&b) distribution can be generated efficiently by just reweighting the leading-order cross 
section. 

The dependency of the x integration limits on the y radiative variable, coupled with 
the definitions of the p and plus distributions, is a significant obstacle to the numerical 
implementation. We may obtain fixed integration limits by a simple change of variables 
x — » x, where x is defined through the relation 



x(x,y) =x(y) +r)(yf x, 



v (y) = V 1 ~ x (y)- 



(4.2) 



Care must be taken in implementing these changes of variables in the NLO differential cross 
section, particularly with regard to the plus and p distributions. When the transformation 



- 12 - 



is complete only plus distributions remain since, unlike x, x extends over the range [0, 1] for 
all y. Numerical implementation of the B (<&b) distribution requires all plus distributions 



be replaced by regular functions; due to the change of variables in Eq.^2| this is now trivial 
since all plus distributions are to be integrated over exactly the domains specified in their 
definition < x < 1. 

After the change of variables, the generation of the iV-body configurations is technically 
carried out in the same way as in Ref. [32]: TV-body configurations are first generated using 
the corresponding leading-order Herwig++ simulation, after which they are reweighted and 



retained with a probability proportional to the integrand of Eq. |4.1| , which is sampled using 
a VEGAS based algorithm [64]. 

4.2 Generation of the hardest emission 

Given an iV-body configuration generated according to B (^b), we proceed to generate the 
largest transverse momentum emission according to the modified Sudakov form factor in 



Eq.|2.3|. The exponent in the modified Sudakov form factor consists of an integral over a 



sum of different contributions, one for each channel, a + b — ► n + c, given by 

Wab &B) = r—— = Hab Cab {x @ ,X e ) , 4.3 

B (9>bJ 2tt x 

where 7i a b is equal to H a b with the plus and p regularization prescriptions omitted. 

Instead of generating the hardest emission in terms of = {x, y} we find it more 



convenient to make a change of variables to & R = {px, y/c}, defined in Eq. 3J3. Making this 



change of variables removes the complicated ^-function in Eq. [2.3| , replacing it by a lower 
bound on the integration over px- 

(/"PTmax \ 
- / d<S> R J2 W ab ($R,$B) , (4.4) 
J PT ab J 

where the upper bound, pr max , is due to the usual basic phase-space considerations. The 
distribution of the transformed radiative variables arising from (pr) (Eq.|2.3|) is sampled 
using a veto algorithm [65], in precisely the same way as was done in Ref. [32]. If an 
emission is generated it is reconstructed from the px and y& radiative variables according 



to Eqs. 3.7,3.5. 



When generating the hardest emission we use a factorization scale of the transverse 
mass of the Higgs boson or off-shell vector boson, in gluon fusion and Higgs-strahlung 
processes, respectively. In both cases we use the px of the boson as renormalization scale. 
This is required to correctly treat the small pt region where the POWHEG results should 
agree with the default Herwig++ parton shower. 

4.3 Truncated and vetoed parton showers 

The Herwig++ shower algorithm [44,66] works by evolving downward in a variable related 
to the angular separation of parton branching products, q, starting at a scale determined 
by the colour flow and kinematics of the underlying hard scattering process, and ending at 
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an infrared cut-off, beneath which further emissions are defined to be unresolvable. Each 
branching is specified by an evolution scale q, a light-cone momentum fraction, z, and an 
azimuthal angle, 4>. The momenta of all particles forming a shower can be uniquely con- 
structed given the {q, z, (ft} parameters of each branching. Since the showers from each 
parton in a given process evolve independently, from what are initially on-shell particles, 
generated according to a matrix element, some reshuffling of these momenta, after the gen- 
eration of the parton showers, is required to ensure global energy-momentum conservation. 

In order to carry out showering of iV+l-body final states associated to the generation 
of the hardest emission, we treat the iV+1 momenta as having arisen from the showering 
of an iV-body configuration with the Herwig++ shower. To this end we calculate the 
branching parameters {q^, Zh, 4>h\ for which the shower would reconstruct the iV+l-body 
system from an initial iV-body one. Details of this inverse reshuffling calculation were given 
already in [32]. The POWHEG emission is subsequently regenerated in the course of a single 
Herwig++ shower as follows: 

1. the shower evolves from the default starting scale to q~h, with the imposition that 
any further emissions be flavour conserving and of lower px than that of the hardest 
emission (pT h ), the truncated shower; 

2. the hardest emission is inserted as a set of branching parameters {q~h, Zh, (fth}', 

3. the evolution continues down to the cut-off scale, vetoing any emissions whose trans- 
verse momentum exceeds pr h , the vetoed shower. 

Should the hardest emission occur in an area of phase space that the shower cannot pop- 



ulate, i.e. the wide angle/high px dead zone (Sect. |5.2. l| ) , subsequent emissions will have 



sufficient resolving power to see the widely separated emitters individually. It follows that 
no truncated shower is then required, since this models coherent, large angle emission from 
more collimated configurations of partons, and so we proceed directly to the vetoed shower. 

5. Results 

In this section we present predictions from our POWHEG simulations of the g + g — > H and 
q + q — > V + H processes within the Herwig++ event generator. By comparing our results 
to other predictions, based on independent calculations and methods, we aim to validate 
these realisations of the formalism. 



In Sect. 5.1 we seek to check the calculation and implementation of the POWHEG NLO 
differential cross section and B(&b) functions, Eqs. [Q] - |2.2| , we thereby check the NLO 
accuracy of the calculation and in particular the generation of the Born variables (Sect. [O]). 
Technically this is the most delicate part of the simulation, requiring a full calculation and 
numerical implementation of the NLO differential cross section. We compare our results to 
the NLO parton level Monte Carlo program MCFM [47] to this end. 

In Sect, |5.2j we move to focus on distributions sensitive to the generation of the hardest 
emission (Sect.^]^) and the subsequent merging with the shower algorithm. Here we com- 
pare our results to three different approaches: the bare angular-ordered parton shower, the 
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parton shower including matrix element corrections and also MCONLO. While MCONLO 
consistently combines NLO calculations with the HERWIG parton shower, matrix element 
corrections work to adjust the distribution of the hardest emission from the Herwig++ par- 
ton shower to be equal to that of the real part of the NLO contribution. Matrix element 
corrections also serve to populate an area of the real emission phase space which the shower 
cannot ordinarily reach, the so-called dead-zone, which we describe fully in Sect. 5.2 . Since 
the effects of the NLO contributions on the normalisation of the results are examined in 



detail in Sect-ISJJ and, moreover, the predictions from the shower, with and without matrix 
element corrections, have leading-order normalisations, in Sect. 5^2 we concentrate on the 
shapes of these distributions. 

Note that we use the notation O (ag) to denote terms of order relative to the leading 
order contribution. 



5.1 POWHEG differential cross sections and B(<&b) 

In order to check the calculation of the POWHEG differential cross section and B(&b) 
functions, Eqs. 24.-2^, total cross sections and differential distributions were compared 



between our POWHEG implementation and the parton level NLO program MCFM [47]. 
Since MCFM computes the effects of fixed order (NLO) corrections to the processes in 
question, these comparisons are performed prior to any showering of the POWHEG NLO 
configurations by Herwig++. 

In carrying out these comparisons both Herwig++ and MCFM were run using the 
MRST2001 NLO [67] parton distribution set via the LHAPDF interface [68]. A fixed, 
constant, factorization and renormalization scale of 100 GeV was used, in order to eliminate 
small variations in the treatment of the running coupling and PDF evolution as a possible 
source of discrepancy. Also, for this part of the validation, we have assumed the Higgs boson 
mass to be 115 GeV. In all cases the total cross sections from MCFM and our POWHEG 
implementation agreed to within 0.5%. 

The gluon fusion process is rather simple in view of the fact that the scalar Higgs 
boson decays isotropically in its rest frame. Consequently the only non-trivial distributions 
to check are the mass and rapidity of the Higgs boson. These distributions are plotted for 
the process 4 gg —* H — > r + r~ in Fig. [j] and everywhere exhibit a high level of agreement 
with MCFM. 

The Higgs-strahlung process is more involved than the gluon fusion process due to 
the intermediate particle having spin-1 and due to our re-using the method in our earlier 
POWHEG work, employing the Kleiss trick [45,46] to generate the NLO corrections indepen- 
dently of the details of the decay of the initially off-shell vector boson. Therefore, to check 
the correctness of this method, one must look closely at the distributions of the final-state 
particles, particularly, on account of propagating spin correlation effects, the vector boson 
and its decay products. 



4 MCFM does not implement the gg — > H — » 77 process. 
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Figure 1: Comparisons of the POWHEG implementation in Herwig++ and the NLO parton level 
code MCFM [47], for the Higgs boson rapidity (yn) and mass (nan) distributions in the process 
g g — > H — > t + t~ . Results on the left are obtained for pp collisions at the Tevatron (y/s = 1.96 TeV) 
while those on the right are for LHC pp collisions (y/s = 14 TeV). 




Figure 2: Mass of the off-shell virtual vector boson in qq — ► HZ and qq — > HW~ processes, in 
the left and right-hand plots respectively. The qq — > i/Z predictions are for 1960 GeV, Tevatron, 
proton-antiproton collisions, while the qq —> HW~ predictions are for 14 TeV, LHC, proton-proton 
collisions. This mass is one of the so-called Born variables in the NLO differential cross section 
(see Sect.lO and Ref. [32]). 
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Figure 3: The polar angle of the electron produced by the vector boson, in its rest frame, in 
Higgs-strahlung processes. The uppermost plots concern the qq — > HW~ process and the lower 
plots relate to the qq — > HZ process. For each process we have displayed the results obtained 
at the Tevatron on the left and the LHC on the right. The leading-order (LO) predictions of 
Herwig++ and MCFM [47] are shown as dashed lines while solid lines represent the corresponding 
NLO predictions. This is a very important test of the correctness of the Kleiss trick, which we have 
used to generate the NLO corrections independently of the generation of the LO process [32]. 



The variety of Higgs-strahlung processes which can be simulated by MCFM is limited, 
unlike Herwig++, hence we opted to carry out the Higgs-strahlung comparisons using the 
following processes: qq — > HW + — > e + v e bb 1 qq — ► HW~ — ► e~v e bb, qq — ► HZ — ► e + e~bb. 
In Fig.^ we show the mass of the initial off-shell vector boson, one of the Born variables 
for this process, for which there is very good agreement between our POWHEG result and 
that of MCFM. In Figures^[]6[ we show a number of distributions with sensitivity to the 
details of the decay of the vector boson, specifically, the polar angle of the lepton produced 
by the decaying, resonant, vector boson in its rest frame, the pseudorapidity of the final- 
state lepton, as well as the rapidity and transverse momentum of the resonant vector boson 
decaying to leptons. In all cases the agreement between our code and MCFM is very good. 
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Figure 4: The pseudorapidity of the electron produced by the decaying vector boson, in 
qq — * HW + (left) and qq — > HZ (right) Higgs-strahlung processes, from the POWHEG imple- 
mentation and also MCFM. As in Fig.|| the predictions for the Tevatron are given on the left 
(pp, y/s = 1.96 TcV) and for the LHC on the right (pp, \/s = 14TeV). As well as the distribution 
of the lepton polar angle, this is a critical test of the implementation of the Kleiss trick which 
factorizes the generation of the LO and NLO calculations [32]. It shows that all spin correlations 
have been correctly propagated through to the vector boson decay products. 



5.2 Hardest emission generation and showering 

In this section we focus on the distributions most sensitive to the generation of the hardest 
emission (Sect. fO| ) and any further radiation obtained due to merging with the shower 
algorithm (Sects. [2], fO| ). In particular we study the pr spectra of the Higgs boson in 
the gluon fusion process and the colourless vector boson plus Higgs boson system in Higgs- 
strahlung processes. The distributions of the rapidity difference between the leading, highest 
Pt, jet and the produced colour neutral systems, yj et — yu and Vj e t — Yvh, are given special 
attention, particularly in view of the differences noted in previous works on the POWHEG 
formalism, arising between it and MC@NLO [28,33,41]. 

Since we generate the hardest emission directly in terms of px it is clear that the px 
spectra are a direct test of this part of the work. The relevance of the yjet— YH and yjet — yvH 
distributions to these investigations is not immediately obvious. However, the yj e t — V H and 
yjet — yvH variables, for one emission, can be expressed purely in terms of the radiative 
variables x and y, they are in fact both equal to y& — y as given by Eq. 3^, hence they 
are also an ideal probe of the hardest emission generation. In order to gain some physical 
insight into what this variable really means, we note that in the limit that the angle between 
the radiated parton and colliding beam partons tends to 90°, in the partonic centre-of-mass 
frame, Eq. [TI^ approximates to 

yk-y = --^- (o-Vi , (5-1) 



1 + x V 2 

where we remind the reader that 9 is the angle between the emitted parton and the p, 
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parton in that frame. Furthermore, from Eq. 3J) one can see that y^ — y is maximised when 
the radiated parton is anticollinear to p§ and minimised when it is collinear to p®. Put 
simply, the central region of the y^ — y distribution is populated by wide angle radiation 
while the tails are due to collinear radiation. 

The jets are defined using the longitutinal invariant algorithm [69] with an angular 
measure AR = \J Ar\ 2 + A<p 2 , where Arj and Acft are the pseudorapidity and azimuthal 
angle differences between two particles respectively, and the E recomination scheme, as 
implemented in the KtJet package [70]. 

In studying the Higgs-strahlung process we shall compare our results to three different 
approaches, namely, the bare angular- ordered parton shower in Herwig++, the Herwig++ 
parton shower including matrix element corrections and also MCONLO. In the case of the 
gluon fusion process we compare to a further fourth prediction which is given by modifying 
the hard component of the matrix element corrections in Herwig++. 



5.2.1 The dead zone 

In order to better understand the predictions from the Herwig++ shower, with and without 
matrix element corrections, and also those of MCONLO, it will help to understand the phase 
space available for the first emission in the shower approach. 

In general, given a leading-order configuration, one must first specify starting scales 
for each parton to evolve down from. The guiding principle behind the choice of starting 
scales in the HERWIG and Herwig++ angular- ordered shower algorithms is as follows: given 
two colour connected shower progenitor partons a and b, progenitor a cannot emit any soft 
radiation into the hemisphere defined by the direction of progenitor b, in the rest frame of a 
and b, and vice-versa. This avoids any double counting of the phase space 5 . In other words, 
the starting scales are fixed by constraining that, in the limit of soft emissions, the two jet 
regions, of the first emission phase space (those which either progenitor can emit into) meet 
smoothly and do not overlap. If one then plots, in the full phase space, the contours to 
which these values of the evolution variables correspond to, one finds three regions: a jet 
region into which progenitor a can emit, a jet region into which progenitor b can emit and a 
further region into which neither can emit, the so-called dead zone. We show exactly these 
phase-space regions for the HERWIG [71,72] and Herwig++ [44,66] algorithms in Fig.[7|, 
taking the case of gluon fusion at LHC energies as an example. 

At this point we wish to remind the reader that in the distributions of yj e t — yH and 
Yjet — yHV which follow, the region around zero corresponds to the emission of the jet at 
right angles to the colliding partons in the partonic centre of mass frame (Eq. |5.1|) . This 
in turn corresponds to the region either side of the line y = in the phase map shown in 
Fig. [?], along which the volume of the dead zone is maximised. 

In Fig.|| we superimpose, on the x, y phase-space map of the gluon fusion process, 
four contours corresponding to constant values of px- 10 GeV, 40 GeV, 80 GeV and mn. 
This was done for the three scenarios which we study using the Monte Carlo predictions, 



6 In practice, in the case of HERWIG, there is a small amount of overlap in the phase space allotted to 
each shower progenitor, although the algorithm later corrects for this by a vetoing procedure. 




Figure 5: The rapidity distributions of the resonant vector boson in the Higgs-strahlung pro- 
cesses, obtained with MCFM and the Herwig++ POWHEG implementation. The left-hand plot was 
obtained considering the qq — > HW~ process while the right-hand plot relates to the qq — > HZ 
process. As with the previous figures, on the left-hand side we show predictions for the Tevatron 
and, on the right, predictions for the LHC. This is another fundamental test of the functioning of 
the Kleiss trick, since, in using this trick, the generation of the radiative variables is completely 
independent of the generation of the decay of the virtual vector boson viz. the Higgs boson and the 
resonant vector boson. 




Figure 6: The transverse momentum of the Z boson in Higgs-strahlung events at the Tevatron, 
compared to MCFM. On the right hand side we show an analogous comparison for the LHC. This 
distribution is a further confirmation that the Kleiss trick is working correctly, since it is sensitive 
to the details of the final state and not simply the production of the initial off-shell vector boson. 



specifically, a Higgs boson mass of 160 GeV at Tevatron energies (1960 GeV), a Higgs boson 
mass of 115 at LHC energies (14 TeV) and finally a Higgs boson mass of 300 GeV, also at 
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Figure 7 : The full radiative phase space is given by the black rectangle in the x, y plane, bounding 
the plots on the right. The minimum kinematically allowed value of x is shown as l/s max = p 2 /s. 
The region into which two colour connected partons a and b cannot emit radiation, according to the 
Herwig++ and HERWIG shower algorithms, is marked 'dead zone' on the left- and right-hand plots 
respectively. In both algorithms, radiation into this dead zone is only possible with the help of a 
hard matrix element correction. Unlike Herwig++, in the case of the HERWIG algorithm there is a 
region of phase space which is double counted by the showering from a and b marked 'overlap'; this 
double counting is ultimately corrected for by a veto procedure. These plots correspond specifically 
to the case of a 115 GeV Higgs boson being produced via the gluon fusion process at LHC energies, 
although their form does not change significantly for the other processes we consider. 

LHC energies. In doing this we see that restricting the phase space to regions with higher 
and higher px leads to the expected result, namely, that the dead zone begins to fill the 
allowed region. 

For the three scenarios we consider, Fig.|| shows that a cut of px >80 GeV already 
leads to a great reduction in the area of phase space populated by the shower, while for 
Pt > rnn the allowed region is fully contained within the dead zone. Similar plots are 
obtained for the Higgs-strahhmg process by simply replacing the Higgs boson mass used in 
the calculation of the gluon fusion phase space, by the typical mass of the colourless vector 
boson plus Higgs boson system (Fig.|2| shows this to be in the range 200-300 GeV). The 
maps in Figs.f?] and |8] are key to a good understanding of the results in Sects. |5~2~2| , |5.2.3 . 

As noted earlier, the dead zone of phase space, into which the shower cannot emit, can 
be filled with the aid of a hard matrix element correction. This involves populating that 
region according to the single real emission matrix element squared [44,73]. In principle 
this is a simple procedure, since the dead zone does not run into any singular regions of 
phase space 6 : given an underlying iV-body configuration one selects whether an emission 
into the dead zone occurs according to the conditional probability, 



and, if an emission is to be generated, it is distributed according to the single real emission 
matrix element squared including PDFs (R (<3?,b, &r))- Neglecting terms beyond NLO ac- 

6 Although the 'throat' of the dead zone in Herwig+- 1- touches the soft boundary, it does so only in a 
vanishing region of the dead zone phase space, which is anyway cut off by the gluon mass regulator used in 
the shower. 




(5.2) 
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Figure 8: Radiative phase space in the x, y plane, as in Fig.^ but with contours corresponding 
to constant values of pr superimposed in green, from lightest to darkest (right to left) respectively 
these are, pr = 10 GeV, pr =40 GeV, pr =80 GeV and pr — vein- The region into which the two 
shower progenitors cannot emit radiation is again marked 'dead zone'; for HERWIG this area is 
bounded in red while for Herwig++ it lies between the magenta lines. The three plots correspond 
to the three scenarios under study in the remainder of the paper: a Higgs boson of mass 160 GeV 
at the Tevatron, a Higgs boson of mass of 115 GeV at the LHC, and a Higgs boson of mass 300 
GeV also at the LHC. The area with emissions pr > ran is entirely within the dead zone in all 
scenarios. 

curacy, V^^ d {&b), integrated over the Born variables, gives the fraction of the NLO cross 
section which the dead zone would contribute. 

In all cases the hard matrix element correction is accompanied by a soft matrix element 
correction, which corrects the distribution of the hardest emission in the parton shower 
regions so that it is also given by R &r). The combination of the soft and hard matrix 
element corrections ensures that, to O(as), the distribution of real radiation is exactly 
matched either side of the dead zone boundaries, i.e. any sensitivity on the position of the 
dead zone boundary will be at the level of O (a|) terms, which should not be logarithmically 
enhanced since the boundary is predominantly in the high pr region. 

The same exact matching at O (as) is true in the MCONLO program, which feeds 
events to the HERWIG parton shower. In this case the generation of the first emission 
is done according to the full NLO differential cross section with additional, resummed, 
higher order corrections entering in the shower regions of the phase space. Sensitivity to 
the dead zone boundary is present in the NLO calculation through the shower subtraction 
terms, which are required to avoid double counting of the NLO contributions through the 
subsequent showering with HERWIG. Ultimately, in MCQNLO, one has that the emission 
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rate in the shower regions is given by the resummed rate in the shower's Sudakov form 
factor, corrected at O (as), while in the dead zone it is related to the fraction which that 
area contributes to the total NLO cross section. This is much like the case of the matrix 
element correction procedure and, as with that method, any mis-match between the shower 
region and the dead zone should therefore again relate to unenhanced O (a|) terms. We 
stress that, unlike the matrix element correction method, the M CO NLO program includes 
also exact NLO virtual corrections to the process in the event generation process, giving 
full NLO accuracy. However, as we now wish to focus on the shapes of the distributions 
sensitive to the effects of real radiation, neglecting the overall normalisation, differences 
arising from virtual corrections only enter the results through terms beyond NLO accuracy. 



This last point is addressed further at the beginning of Sect. |5,2,2 , 

As we have already described in Sects. |2| and ||, the POWHEG method generates the 
hardest emission completely independently of the detailed workings of the shower, it may 
be considered as being a pr ordered shower in its own right, albeit for a single emission. 
Unlike the other methods, POWHEG therefore has no dependence whatsoever on the dead 
zone boundary and the way in which it generates the hardest emission generation is the 
same throughout the whole phase space. 

5.2.2 Gluon fusion 

In this subsection we compare our POWHEG simulation against predictions from the bare 
angular-ordered parton shower in Herwig++, the Herwig++ parton shower including matrix 
element corrections and MC@NLO. We also compare against a further prediction from the 
matrix element correction (MEC) procedure, in which we decrease the rate of emission into 
the dead zone by unenhanced terms of O (a|) , by changing the denominator in V^™ d (3>b) 
from B($ B ) to 

« (*b) - «d° (*b) = / d* R • (5-3) 

./dead -D (3>£ J 

Whereas integrating P^ead (^b) over ^ ne Born variables gives the fraction of the NLO cross 
section contributed by the dead zone neglecting higher order terms, performing the same 
integral with pjg£ ($ B ) gives this fraction exactly 7 . Although it should be obvious to the 
reader that this is technically only an alteration at the level of terms beyond NLO accuracy, 
we do not wish to give the impression that it is a priori a small change, at least not for 
the gluon fusion process (recall Fig.[j]). We expect that this alteration should mean that 
the modified MEC predictions should reproduce well the rate at which MC@NLO emits 
radiation into the dead zone. 

In figure ^ we show the px spectrum of the Higgs boson and also that of the hardest jet. 
One can see that the pr spectra at the Tevatron are less hard than those obtained at the 
LHC, as one would expect given the greater centre-of-mass energies of the latter. One also 
expects that, at the LHC, the larger mass of a 300 GeV Higgs boson would automatically 
give rise to it having a harder spectrum than that of a 115 GeV Higgs boson, this is indeed 
the case for all five Monte Carlo predictions. 

7 For this reason we choose to distinguish the modified probability by the superscript NLO. 
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Figure 9: Transverse momentum spectra for the Higgs boson and the leading, highest pr, jet, 
obtained using Herwig++ with matrix element corrections (black), Herwig++ without matrix ele- 
ment corrections, i.e. the uncorrected parton shower (dotted), MCONLO (blue) and our POWHEG 
simulation inside Herwig++ (red). The green curve is obtained by modifying the hard component 
of the matrix element corrections, to decrease the amount of radiation produced in the associated 
high pr, wide-angle, dead-zone, by terms beyond next-to- leading-order accuracy. This modification 
is discussed further in Sect. 5.2. 



- 24 - 



In each of our three phenomenological scenarios we see that the behaviour of the dif- 
ferent simulations is more-or-less the same with respect to one another. The effect of the 
radiation dead zone is very clear as a sharp knee in the spectra from the uncorrected Her- 
wig++ parton shower (black dotted lines), as the transverse momentum approaches the 
Higgs boson mass. Only the effects of multiple emission mean that this is not an abrupt 
cut off (Fig.|8|). As predicted, turning on the MECs and hence filling the dead zone in 
Herwig++, leads to a spectrum in much better agreement with all of the other methods 
(black lines). A very slight kink is still visible in this default MEC prediction, which we 
already forecast in Sect. 5.2.1 as a mis-match, at O (a|), across the dead zone boundary. 



The modified MEC (green dotted lines) shows a similar trend with respect to the 
uncorrected parton shower prediction, although it comes as no surprise that it has a softer 
spectrum than the regular MEC, simply because the rate of emission into the dead zone is 
reduced by an amount approximately given by the NLO K-factor cf. Eqs.^1, 5.3. From 



our earlier investigations concerning the dead zone one should expect that the difference in 
the emission rates in Eqs. 5^ and 5J3 directly manifests itself as the same relative difference 
in the upper values of the px spectra of the Higgs boson and the leading jet. This is indeed 
seen to be the case in Fig. [9], where the unmodified MEC prediction is around a factor of 
two higher than the modified one on the right-hand side of each plot. 



In all cases the Herwig++ NLO POWHEG prediction (red lines) is seen to be in very 
close agreement with the predictions obtained using the normal MEC procedure (black 
lines). As we shall discuss more in Sect. 5^5, we attribute this to the fact that the emission 
rates for the MEC method and the POWHEG implementation can be expected to converge 
to the same value at high px- 



The M CO NLO px spectra are in general somewhat softer than those of the MEC pro- 
cedure and POWHEG NLO prediction, they tend to lie between the prediction obtained 
using the modified MEC (green dotted lines) and the Herwig++ POWHEG/default MEC 
lines. Good agreement between the MC@NLO and modified MEC predictions can be seen 
for the case of a 115 GeV Higgs boson at the LHC, similarly, the agreement in the case 
of the Tevatron is good below 200 GeV, however, for a 300 GeV Higgs boson at the LHC 
the agreement is not as good as hoped beyond px — w-h- We suggest that these discrep- 
ancies arise from presence of an improvement introduced between versions 3.2 and 3.3 of 
MC@NLO, with the aim of improving the description of the Higgs boson px spectrum [22]. 



In figures 1C and [Tl] we show the distribution of the rapidity difference between the 
leading jet and the Higgs boson, for increasingly hard cuts on the px of the jet. Starting 
again with the bare parton shower prediction (black dotted lines), for each scenario and each 
value of the px cut on the leading jet, the structure is broadly the same, the distribution 
rises from the tails at either side of the plot into a hump before falling again into a very 
deep dip in the centre, at yj e t — yH = 0. 
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Figure 10: In the left hand column we show distributions of the difference in rapidity between 
the leading jet and a Higgs boson of mass 160 GeV at the Tevatron, for increasing cuts on the pr 
of the leading jet. For a single emission, in the region close to zero, this variable is proportional to 
the angle between the emitted parton and the transverse direction in the partonic centre-of-mass 
frame (Eq. 5A). In the right hand column we show the corresponding jet multiplicity distributions. 
The black and dotted lines show the predictions obtained using Herwig++ with and without matrix 
element corrections, respectively. The blue line shows the prediction from MC@NLO and the red 
line is that of our POWHEG simulation in Herwig++. The green line is obtained using a modified 
version of the hard matrix element correction, effectively decreasing the amount of radiation that 
this method produces in the high px, dead zone, by terms beyond NLO accuracy (see Sect. 5.3). 
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Figure 11: Distributions of the difference in rapidity between the leading jet and the Higgs boson 
in the gluon fusion process at the LHC, for increasing cuts on the pr of the leading jet. The series of 
plots on the left hand side are obtained for a Higgs boson of mass 115 GeV, while those on the right 
correspond to a Higgs boson of mass 300 GeV. The colour assignment of the various predictions is 
described inset, it is the same as for earlier Tevatron predictions in Figjtol 
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These volcano formations are a direct manifestation of the radiation dead zone; the 
volume of the dead zone is maximised along the direction corresponding to wide angle 
radiation in the partonic centre-of-mass frame, y « 0, 8 « | (Sect. 5.2.1 , Figs.[?], |8|) and the 
Yjet - Yh variable is proportional to 9 — \ (Sect.[T^, Eq. |5.1| ). The only thing which may 
then be slightly surprising is that the distribution does not in fact go exactly to zero at 
yjet — Yh = 0, this is simply due to the effects of multiple emissions (the phase space in 
Figs.^ and |S] is only exact for the case of a single emission). 

Looking in more detail at the structure of the uncorrected Herwig++ predictions, we 
see that the throat in each of volcano distributions widens as the pt cut on the hardest 
jet is increased. This broadening occurs because the dead zone accounts for proportionally 
more and more of the smaller/less wide angle regions of the phase space as the pt increases. 
This trend of the broadening throat with the hardening of the pt cut is therefore visible in 
each of the three scenarios we consider and it will come as no surprise that the same occurs 
for the Higgs-strahlung process. 

The hard matrix element correction begins to fill in the central throat region by emitting 
wide angle radiation into the dead zone. Nevertheless, both the modified and unmodified 
MEC results carry some residual sensitivity to the dead zone boundary, as evidenced by 
their inheritance of dips in the central region, dips which follow the same trend as that set 
by the uncorrected parton shower, to broaden and deepen as the pr cut on the hardest jet 
is increased. We attribute this dip behaviour as being almost entirely due to the mis-match 
in the emission rate across the dead zone boundary, which occurs for terms beyond NLO 
accuracy, as discussed in Sect. 5.2.1 . This conclusion may not appear to hold for the 115 
GeV Higgs boson when the pt cut on the jet reaches 80 GeV. However, in this plot, the 
absence of a dip actually reinforces our assertion, as can be seen by consulting the 80 GeV 
Pt contour in the corresponding phase-space map; this contour shows that the allowed 
region for emissions is almost identically in the dead zone, so any matching across the dead 
zone boundary and, conversely, any mis-matching, is extremely limited in this special case. 

Having understood how the dead zone is manifest in the rapidity difference plots, 
the differences observed between the modified and unmodified MEC methods are rather 
unremarkable: since the emission rate of the modified MEC into the dead zone is reduced 
by around a factor of two with respect to the unmodified correction, the former emits more 
wide angle radiation and so populates the central region of the yj et — yH to a greater extent. 

Although the HERWIG program has a slightly different phase-space coverage for its 
parton showers (Figs. [7], ||), they are apparently not so different, particularly away from the 
soft region. Considering the region of the phase space allowed by just a 10 GeV pt cut 
(Fig.||) it is already clear that the coverage by HERWIG and Herwig++ is really very similar 
and that it is basically identical by the time the cut reaches 80 GeV. The predictions of the 
uncorrected HERWIG shower, which showers the positive and negative unit weight events 
fed to it from MC@NLO, will then be very similar to those shown here for Herwig++, in 
particular, the volcano structures in the yj e t — yH distributions. This being the case, it is 
understandable that the M CO NLO distributions (blue lines) also exhibit dips in the central 
region, and that the behaviour of these dips with the varying px cut follows that of the 
MEC method; we also note that a number of these results are markedly similar to those 



- 28 - 



obtained using the modified MEC. 

Before discussing the POWHEG results we reiterate that this method is wholly indepen- 
dent of the details of the partitioning of phase space in the shower algorithms, it generates 
the NLO emission effectively as a self contained pt ordered parton shower, albeit with NLO 
accuracy, and thereby circumvents such issues. Hence, the appearance of a slight dip in 
the distributions of Vj e t — yH cannot be explained in the same way as those seen in the 
predictions of the Herwig++ shower, with or without the MECs, nor those of MCONLO. 
The central dip appearing in the POWHEG results alters the height of these distributions, 
in all cases, by less than 5%, it is therefore characteristic of the uncertainties typical of 
NNLO calculations, and is many times smaller than those seen in the other predictions. 
Given the smallness of the effect, its absence in the Tevatron distributions, and also the 
fact that it exhibits no discernible response to the changing pt cut, we cannot comment on 
its origin. 

Finally we note that only the POWHEG predictions clearly exhibit the expected physical 
behaviour on increasing the pt cut: that the yj e t — yH distribution should become squeezed 
toward the central region, as the phase space available for small angle emissions, which 
populate the tails, becomes reduced relative to the phase space available for wide angle 
emissions. This trend is somewhat obscured in the distributions predicted by the other 
methods. 

At this point we do not wish to give the impression that the dips exhibited by the MEC 
and M CO NLO predictions are incorrect. It is our contention, however, that these predictions 
can be consistently explained by ascribing the origin of the dips to a mis-match at O (a|) in 
the emission rates either side of the dead zone boundary. Although the dead zone boundary 
is an unphysical partition in the phase space, we stress that the mis-match involves terms 
beyond the stated accuracy of either method. Looking ahead, beyond NLO accuracy, it 
should then be the case that the MEC and M CO NLO approaches fail to approximate any 
higher order terms, while the POWHEG method, being independent of any kind of artificial 
phase space partitioning, should fare much better. 

In figure [10| we also show jet multiplicity distributions associated to each of the Tevatron 
Yjet — Yh plots. As expected, these plots show, in all cases, that the jet multiplicity distribu- 
tion decreases rapidly as the px cut on the leading jet is increased, with only ~5% of events 
containing a jet once the cut reaches 80 GeV. In the case of the soft pt cut (pr>10 GeV) 
one can see that the POWHEG approach predicts events with lower multiplicity than the 
other methods (which broadly agree with one another). 

5.2.3 Higgs-strahlung 

In figure [l^ we compare the transverse momentum spectra of the Higgs boson and W 
boson assuming a 160 GeV Higgs boson mass at the Tevatron and Higgs boson masses of 
115 GeV and 300 GeV at the LHC. In each case we compare the results obtained using the 
uncorrected Herwig++ parton shower, the parton shower with MECs, MCONLO and our 
POWHEG implementation. All four approaches agree remarkably well. The fact that the 
uncorrected, leading-order, parton shower prediction agrees so well with the other methods, 
which include at least the NLO real emission corrections, indicates that these distributions 
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are rather insensitive to the emission of additional radiation, therefore one should not expect 
to see differences among the more sophisticated approaches. 



On the right of figure 13 we show the rapidity of the leading jet in qq — > ZH events. 
These rapidity distributions are concentrated more in the central region in the case of the 
Tevatron than the LHC. This behaviour can be inferred from the fact that there is more 
phase space available for extra radiation at LHC energies. The same line of reasoning also 
explains why the rapidity distribution in the case of the 300 GeV Higgs boson at the LHC 
is also more central than that of the 115 GeV Higgs boson. There is a tendency in all 
of the plots for the MCONLO distributions to contain more events in the tails, conversely, 
the POWHEG results show more jets produced in the central region. The predictions from 
the uncorrected parton shower and the parton shower with a MEC lie between those of 
MC@NLO and POWHEG, with the former being slightly closer to MCONLO and the latter 
closer to POWHEG. 



The plots on the left of figure are a more interesting test of our methods, they show 
the transverse momentum of the colourless Z-Higgs boson system and also the rapidity of 
the hardest jet in the same qq —> ZH events. The pr of the vector boson plus Higgs boson 
system is generated directly by the Higgs-strahlung POWHEG simulation. As in the case of 
the gluon fusion process we see that the POWHEG results and those of the parton shower 
including matrix element corrections are essentially the same, and that both are harder than 
the corresponding MCONLO prediction. However, the degree by which the latter predictions 
are above those of MCONLO is significantly reduced with respect to that seen in the gluon 
fusion case. We again attribute this to the relative differences in the rate of emission into 
the dead zone, in the MEC method this occurs with probability Vf^ d (Eq.fJ) while 
in MCONLO the analogous probability will be essentially given by the fraction which the 
dead zone contributes to the total NLO cross section i.e. (Eq.pD. Since the 

denominators in these emission probabilities differ by an amount of order the NLO K- 
factor the differences arising in the gluon fusion case should be large whereas they should 
be small in the Higgs-strahlung case. It is difficult to see it clearly at the upper end of 
these Higgs-strahlung pr spectra but the MCQNLO result is below that of POWHEG and 
the MEC methods by roughly 30%, which is compatible with the enhancement due to the 
NLO corrections seen in the comparisons with MCFM in e.g. Fig.^. 



In figures 14 and ^ we show distributions of the rapidity differences yj e t — yzH and 
Yjet — Ywh, for the qq — > WH and qq — > ZH processes, at the Tevatron and LHC respec- 
tively. These plots display the same features and trends as seen in the gluon fusion case. As 
before, the uncorrected parton shower predictions give rise to volcano shaped distributions 
due to the dead zone in the phase space and increasing the pt cut on the leading jet again 
has the effect of broadening the throat, reflecting the fact that the dead zone occupies an 
increasingly large fraction of the allowed phase space. 
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Figure 12: Transverse momentum spectra for the Higgs boson (left) and the W ± boson in 
qq — > HW ± events obtained using Herwig++ with and without matrix element corrections (black 
and black dots respectively), MCONLO (blue) and our POWHEG simulation inside Herwig++ (red). 
The first uppermost predictions, for the Tevatron, are obtained assuming a Higgs boson mass of 
160 GeV. The following four plots are analogous projections for the LHC for Higgs boson masses of 
115 GeV and 300 GeV. 
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Figure 13: On the left we compare predictions for the transverse momentum spectrum of the 
colourless final-state system comprised of the Z and Higgs boson (left) in qq — > ZH events. From 
top to bottom, respectively, these results are obtained for Tevatron energies with a Higgs boson 
mass of 160 GeV, LHC energies with a Higgs boson mass of 115 GeV and LHC energies with a 
Higgs boson mass of 300 GeV. On the right hand side we show the corresponding distributions for 
the rapidity of the leading, highest pr, jet. The colour assignment for the different approaches used 
is shown inset, and is the same as that in Fig.O. 
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The MEC and MCONLO methods emit radiation into the dead zone, filling the throat 
region in the yj e t — Ywh and yj e t — Yzh plots, however, they still exhibit a clear sensitivity 
to the dead zone boundary which is manifest as irregularities in the central regions of the 
distributions. For the case that the pt cut on the leading jet is soft, pt >10 GeV, the MEC 
predictions show a tiny, sharp dip around the centre, while those of MC@NLO show the 
formation of a small tower in the same place. This is plainly a mis-match of O (a|) terms 
across the phase-space partition. As the px veto increases the dips in the MEC predictions 
tend to deepen and in the case of MCQNLO, the small towers turn to small dips. In all 
cases the POWHEG distributions are smooth, exhibiting no irregular volcanic features, as 
expected, furthermore they are more concentrated in the central region than all of the other 
predictions, indicating a tendency to emit proportionally more wide angle radiation. 

We point out that in all of the MEC and MCONLO predictions the residual effects of 
the phase-space dead zone are felt much less strongly in the case of Higgs-strahlung than 
in gluon fusion. These behaviours are less marked for the same reason that the MCONLO 
and POWHEG pr spectra agree better for Higgs-strahlung than for gluon fusion: the fact 
that the NLO corrections are substantially smaller for Higgs-strahlung means B($>b) is 
less different to B ($b) and so the O differences in the rates at which each method 



populates the dead zone are greatly reduced, cf. Eqs.5.2, 5.3 



5.3 Emission rates in the dead zone 



In this subsection we present an heuristic discussion of the rates with which each approach 
emits into the high pt region. In particular we consider the area of phase space correspond- 
ing to transverse momenta pr > m n , where m n is the mass of the colourless final-state 
system. This region is completely contained within the dead zone and the contribution 
which it makes to the cross section is not logarithmically enhanced (Sect. [5~2?i" ). 



From Eq. it follows that the probability for Herwig++ to generate an emission with 
Pt > m n is 



VT (*b) 



d<5> Rl 



B($ 



B 



(5.4) 



where the omitted higher-order terms in the exponential series are negligible, not simply 
because they carry higher powers of the coupling constant but also, significantly, because 
the contribution to the total cross section from this region is in general very small. In Eq. 5.4 
R± ($£?, ) is, as before, the real single-emission matrix element squared, including flux 
factors and PDFs, and and are the Born and radiative variables parametrising 
the two-body phase space. The lower limit on the integral signifies that it extends over 
all available phase space above px > rn n . Our adding a subscript '1' to each R, has no 
significance here but it will be useful later, in discussing NNLO corrections. 
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Figure 14: In the left-hand column are distributions of the difference in rapidity between the 
leading jet and the colourless W— Higgs boson system, in qq — > HW ± events at the Tevatron, 
for increasing cuts on the pr of the jet. For a single emission, the central value yj et — vhw = 
corresponds to a configuration where the 14^-Higgs boson and the colliding partons travel at right- 

In the plots on the right-hand side 
The black and dotted lines show the 



angles in the partonic centre-of-mass frame (Eqs. |3,9| , 
we show the corresponding jet multiplicity distributions, 
predictions obtained using Herwig++ with and without matrix element corrections respectively. 
The blue line shows the predictions of MC@NLO and the red line corresponds to our POWHEG 
simulation inside Herwig++. 



- 34 - 




Hardest jet rapidity - HZ rapidity (p T >40 GeV) Hardest jet rapidity - HZ rapidity (p T >40 GeV) 




Yjet YhZ Yjet YhZ 




Figure 15: Distributions of the difference in rapidity between the leading jet and the Higgs boson 
in the qq — > ZH process at the LHC, for increasing cuts on the pt of the leading jet. The series of 
plots on the left hand side are obtained for a Higgs boson of mass 115 GeV, while those on the right 
correspond to a Higgs boson of mass 300 GeV. The colour assignment of the various predictions is 
described inset, it is the same as for earlier Tevatron predictions for qq — > HW^ in FigjlJ. 



- 35 - 



The fraction of MC@NLO events with emissions in this region is given by the corre- 
sponding fraction of the NLO cross section 

= / d* fll gli^lgl) . (5.5) 



For what follows it will be useful to note that if we expand the denominator in Eq.pJj 
neglecting terms O (a|) and above, we have 



,NLO 



Vm - { * b) = L Ri b{* B ) [ 2 " bWs)) ■ (5 - 6) 



In POWHEG the Sudakov form factor, A^ (px) in Eq^, is the probability that no 
radiation is emitted from the leading-order partons in evolving from the maximum kine- 
matically allowed transverse momentum, down to px- Hence a POWHEG simulation will 
emit radiation into the same high px region with probability V™ = 1 — A^ (m n ). Neglect- 
ing kinematically suppressed 0(a|) terms, as in the case of the MEC, this gives, 

Vm *-L d ** b(* b) • (5 - 7) 



Comparing Eqs.^J and p.Tj we see that the emission rates in the high px region, from 
POWHEG and the MEC, are the same up to the scale used in the strong coupling constant 
(in the former it is px while in the latter it is my). These corresponding emission rates 



are evident in the px spectra in Figs.H and 13. It is also clear that MCONLO will emit 



into this region less often than the other two since, whereas the POWHEG and Herwig++ 
simulations have an emission probability inversely proportional to the Born cross section 
.B (<&#), the probability of emission in MC@NLO is inversely proportional to the, larger, 
NLO cross section, B ($b)> This argument is also supported by what we have seen in the 
px spectra. Essentially we have 

V™*V™*KVZ° (5-8) 

where K, is the relevant NLO K-factor. For processes with smaller K-factors the differences 
in the population of the dead zone in each approach should not differ too much but M CO NLO 
will be systematically softer than POWHEG and Herwig++. 

Is one of these probabilities more correct than the others? In Ref. [41] the authors 
compared the Higgs boson px spectrum from their POWHEG simulation to NNLO predic- 
tions [56,57] and found that they agreed better than those of MC@NLO. The amount of 
radiation produced in the high px dead zone in the case of the fixed order, NNLO, Monte 
Carlo used in Ref. [41] should be the relevant fraction of the NNLO cross section 

J m n L J 

+ dCTNNLO , (5-9) 

where the sum of the last two terms in the numerator represents the finite combination of 
the one-loop single emission matrix element interfering with the tree-level single emission 
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matrix element, and the squared double emission matrix element respectively, including 
PDF and flux factors. Since the numerator of V^ LO is O (as), if we omit terms O (a|) in 
^m!! LO we can replace the NNLO differential cross section in the denominator by the NLO 
one (Eq. 3.24| ) , whence it follows that 



V 



J m„. 



Rl (* B) *Jii) 



1 B($ B ) | RtjOfB^Ri, 



B($b) R x {* Bi $Ri) 



(5.10) 



where B (&b) was defined in Eq.[2^ and i? ^Ri) is defined analogously as 

Si(*B,*fl x ) = ^l($B,^fl 1 )+-Rl+l($B,^ 1 )+ / d^lfe^B,*^!,*^) -(5.11) 



Whereas the NNLO rate has R\ in the numerator all of the others just contain the 
term corresponding to the squared single emission matrix element R\. Replacing i?i — > i?i 
in Eq. 5.10| gives Eq. 5ij. We can consider B/B and R\/R\ as K-factors for the processes 
gg —* i7 and 55 — ► -ff+jet, differential in <I>b and {^b, ^Ri}, respectively. One can see from 



Eqs.lSJ] and p,10| that if these two K-factors coincide so too will the rates for the POWHEG 
and NNLO calculations. In fact this seems to be more-or-less the case for gg —* H and 
gg — > H + jet processes, which have very similar K-factors of around 1.6 / 1.7 [74]. It then 
seems quite feasible that the pt spectrum of the Higgs boson in gluon fusion should be 
better modelled by the POWHEG approach, as was found to be the case in Ref. [41]. 

Equation 5.10| does not tell us that POWHEG will generally reproduce higher orders 
better than M CO NLO. According to Eq. 5. 10| if R\/Ri is more-or-less one and B/B is 



significantly greater than one then MCONLO will give a better estimate of "P^ LO than the 
POWHEG approach. 



6. Conclusion 

In this work we have fully realized the POWHEG NLO matching prescription for Higgs 
boson production via gluon fusion and Higgs-strahlung processes, within the Herwig+- 1- 
Monte Carlo event generator, including truncated shower effects to correctly include colour 
coherence. 

The cross sections and parton level NLO distributions were found to be in very good 
agreement with the MCFM NLO Monte Carlo. The shapes of the emission spectra from 
the full simulation, including parton shower effects, are seen to broadly agree well with the 
older matrix element correction method and also M CO NLO. 

We observe that the px spectra from the M CO NLO program tend to be softer than 
those of POWHEG and the matrix element correction method. We ascribe this effect to 
differences, at the level of O (q|) terms, beyond the stated accuracy of all approaches, in 
the rate at which they emit radiation into the so-called dead zone, associated to high pr 
and wide angle emissions. We have been able to estimate the magnitude of the relative 
difference with a fair degree of success and, based on this line of reasoning, we presented 
an argument for why the POWHEG method appears to better reproduce the NNLO Higgs 
boson pt spectrum in gluon fusion. 



-37- 



Separately, we have shown that a mis-match between terms of this magnitude manifests 
itself as marked sensitivity to the unphysical dead zone partition, for the predictions of the 
rapidity difference distributions yj e t — yH and yj e t — yvH, for both the matrix element 
correction method and MCONLO. These distributions acquire irregularities in the central 
region corresponding to wide angle emissions. This would seem to rule out the possibility 
of these two methods giving a good approximation of unenhanced NNLO effects in these 
distributions. Conversely, the POWHEG predictions for these distributions are smooth and 
physical in appearance, by construction, they have no sensitivity to the dead zone partition. 

At present there is no known method of consistently including exact NNLO fixed order 
calculations within a parton shower simulation. An approximate means of doing this, 
involving the reweighting of PYTHIA [65] and MC@NLO events such that they reproduce 
certain NNLO distributions, has been recently carried out in Ref. [75]. Based on the 
findings discussed above, it is our expectation that applying the same reweighting technique 
to POWHEG events should lead to further improvements in those predictions. 

Both the gluon fusion and Higgs-strahlung simulations we have presented here are al- 
ready available in Herwig++2.3. The algorithm we have used to implement the POWHEG 
formalism, specifically that part concerning the inverse mapping of the hardest emission 
kinematics to a set of shower variables has further applications in multi-leg matching pro- 
cedures e.g. the CKKW method [35,36]. Likewise, the general formulae leading to Eqs.[ 



3.24, for the phase space and differential cross-section for a + b — > n, can be readily used 



to implement other NLO calculations such as those in Refs. [59-61]. This work is already 
at an advanced stage and will appear in the forthcoming version of Herwig++. 
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A. Regularized and unregularized splitting functions 

This appendix contains the splitting functions needed in the NLO calculations. The bare 
(spin averaged) splitting functions in n = 4 — 2e dimensions are of the form 

P 7- (x; e) = P. ~ (x) + eP.V (x) , 

1,IC ' ' 1,IC \ ' i t ic y ' ' 

where 
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— h + X (1 - x) 

1 — X X 



Pgg (X) = 2C A 

P qq (x) = C F 
P qg (x) = C F 

Pgq (x)=T R [l-2x(l-x)], 



1 + x 2 
1 — X 



i + (i-xy 



P e gg (x) = 0, 

P* q (x) = -C F (l-x), 

Pqg ( x ) = ~C F X, 

P* gq (x) = -2T R x(l-x) + 0(e) . 



We write the 'customary' regularized Altarelli-Parisi equations kernels using the ^-distributions 



as 

where 
PP 9 (x) = 2C A 

pP q (x) = C F 

P q P g (X) = C F 



Pile (*) = P-T, 00 + £ [V llc + 41nr/U (1 - x 



x 1 — X 

H h x (1 — x) 



(l-^ 
1 + x 2 



C 99 — Ca, 

Cqq = C F , 



2vr6 
~Ca~ 



3 

^9<? — 2' 



pP q {x)=T R x 2 + {l-x) 



and 

with all other p. ~ and C ~ being equal to zero. 

' ' 1,1c 1,1c c> n 
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